function [log_likelihood, nan_pct]=ml_risk_score_types_cdf(theta)
%%% calculates the log-likelihood allowing multiple factors. Takes
%%% simulated factors from step0 as input 

%%% Works for CARA coeffs without modification (so long as the thresholds,
%%% scale of sig_th, .. are properly adjusted in the master file)


 global X_ind_s scores ind_id factor_n  ind_n draw_n X_n theta_xy lottery_n lottery_choice_ind n_type sim_factor_contribs maxima minima sig_th_scale sig_th_lower_limit 
    

    
    %%% simulated factors from step0
    sim_factors=scores;
    
    %%% Probability contribution from previous step (likelihood of measures
    %%% observed for each individual given his characteristics and
    %%% simulated draws)
    sim_l_contrib_all_factors=sim_factor_contribs;
    
    % # of observations including rows with simulated factor draws for each
    % person
    rows=size(sim_factors,1);
    
    sex=X_ind_s(:,1);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Lottery Part %%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% All Types  %%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%% parameter assignment for risk analysis
    
    theta1=theta(1);
    theta2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];        

    kappa1=theta(1);
    kappa2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];     

    sig_th1=theta(1);
    sig_th2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];
    
  
    type1_3=NaN(3,n_type);
    p_type=zeros(n_type,1);
    
    sim_l_contrib_lottery=NaN(ind_n*draw_n,lottery_n);
    av_l_part=NaN(ind_n,n_type);
    theta_i=NaN(ind_n*draw_n,n_type);
    kappa_i=NaN(ind_n*draw_n,n_type);
    kappa_i_OLS=NaN(ind_n*draw_n,n_type);
    
    %%% weigh over unobserved types, each has its own set of intercepts for
    %%% the structural parameters of interest
    for tt=1:n_type
        %%% last type's proba is just 1-sum(other probas) so don't calculate
        if tt<n_type
            type1_3(:,tt)=theta(1:3);
            %%% probability needs to be between 0 and 1
            p_type(tt)=normcdf(theta(4));
            theta(1:4)=[];      
        elseif tt==n_type
            type1_3(:,tt)=theta(1:3);
            theta(1:3)=[];  
            p_type(tt)=1-sum(p_type);
        end
        
        %%% individual risk aversion based on sex and particular factor draws -
        %%% deterministic part - note that as above, the factor is determined
        %%% by the orthogonal draw AND by X*alpha - a function of observables

        %%% Even though the following parameters are individual-specific - thus
        %%% _i - they vary here within individuals due to the different factor
        %%% draws used in each row for the factor simulation
        theta_i(:,tt)=type1_3(1,tt)+sex*theta1+sim_factors*theta2_n;
        
        theta_i=min(maxima(1), theta_i);
        theta_i=max(minima(1), theta_i);

        %%% individual trembling based on sex and particular factor draws -
        %%% 0.5*normcdf() to bind it between 0 and 0.5 now (one can
        %%% make mistakes at most 50% of time, otherwise, probably
        %%% preferences misestimated and kappa takes on their role).
        kappa_i(:,tt)=0.5*normcdf(type1_3(2,tt)+sex*kappa1+sim_factors*kappa2_n);
        kappa_i=min(maxima(2), kappa_i);
        kappa_i=max(minima(2), kappa_i);
        
        %%% individual variability of risk aversion based on sex and particular
        %%% factor draws 
        sig_th_i=sig_th_scale*normcdf(type1_3(3,tt)+sex*sig_th1+sim_factors*sig_th2_n)+sig_th_lower_limit;
        sig_th_i=min(maxima(3), sig_th_i);
        sig_th_i=max(minima(3), sig_th_i);
        
        %%% compares an individual's level of risk aversion given each set
        %%% of simulated factor draws to the threshold level of risk
        %%% aversion proper to a given lottery choice task
        theta_xy_comp=kron(ones(ind_n*draw_n,1), theta_xy);
        theta_diff=(bsxfun(@plus,theta_xy_comp,-theta_i(:,tt)));
        temp=zeros(size(theta_diff));
        for i=1:lottery_n
            temp(:,i)=theta_diff(:,i)./sig_th_i;
        end
        %%% Probability of choosing the risky option i.e. probability that the
        %%% individual's risk aversion coeff is below the threshold for a
        %%% particular lottery. The difference in the coeffs is above
        %%% divided by the individual's sd of the error term so that the
        %%% normal cdf can be used here.
        %%% assumes normal errors on the risk aversion parameter:
        %%% theta~N(0,(sig_th)^2)
        p_risky=normcdf(temp);
        p_safe=1-normcdf(temp);

        % kappa percent of the time and individual makes the wrong choice.
        choice_risky=zeros(size(p_risky));
        for i=1:lottery_n
            choice_risky(:,i)=p_risky(:,i).*(1-kappa_i(:,tt))+p_safe(:,i).*kappa_i(:,tt);
        end
        choice_safe=1-choice_risky;
        sim_l_contrib_lottery=(choice_risky.^(lottery_choice_ind)).*(choice_safe.^(1-lottery_choice_ind));
        
        
        
       %%%%%%%%%%%%%%%%%%%%%% Combining the likelihood contributions from the various
       %%%%%%%%%%%%%%%%%%%%%% types of measures and lottery task choices
       sim_l_contrib_total=[sim_l_contrib_all_factors sim_l_contrib_lottery];
       % for numerical estimation purposes
       sim_l_contrib_total_plus=sim_l_contrib_total+1*10^(-20);

        % Individual likelihood contribution conditional on factor draws
        sim_ind_draw_contrib=prod(sim_l_contrib_total_plus,2);
        % individual likelihood of measures and lottery choices averaged
        % over simulated factor draws
        av_l_contrib=cell2mat(accumarray(ind_id,sim_ind_draw_contrib,[],@(x){sum(x)}))/draw_n;

        % In order to avoid zeros for the log
        av_l_part(:,tt)=av_l_contrib+1*10^(-200); 
        
        sim_l_contrib_lottery=NaN(ind_n,lottery_n);
        clear sim_l_contrib_total sim_l_contrib_total_plus sim_ind_draw_contrib av_l_contrib sig_th_i theta_xy_comp theta_diff p_risky p_safe choice_risky choice_safe

    end

    %%%%%%%%%%%%%%%%%%%%%%% Combining the likelihood contributions of the
    %%%%%%%%%%%%%%%%%%%%%%% various unobserved types and weigthing them by each type's
    %%%%%%%%%%%%%%%%%%%%%%% estimated population prevalence
    
    av_l=zeros(ind_n,1);
    for tt=1:n_type
        av_l=av_l+p_type(tt)*av_l_part(:,tt);
    end
    
    log_likelihood=-nanmean(log(av_l));
    
    %%% will capture the % of nan and inf in the final average likelihoods
    nan_pct=1-mean(isfinite(av_l));
    if nan_pct>.05
        log_likelihood=NaN;
    end
    
    %%% If at end the calculated proba of the last type is still <0,
    %%% discard estimate
    p_last=p_type(end);
    if p_last<0
        log_likelihood=NaN;
    end
    

    
end